home *** CD-ROM | disk | FTP | other *** search
/ Shareware Grab Bag / Shareware Grab Bag.iso / 090 / byt1186b.arc / MATINV.BAS < prev    next >
BASIC Source File  |  1980-01-01  |  3KB  |  67 lines
  1. 10 REM Real-entry matrix inversion (Pan-Reif method)
  2. 20 REM Program by Thomas E. Phipps, Jr.
  3. 30 REM Ref= Proc. 17th Annual ACM Sympos. on Theory
  4. 40 REM of Computing, Providence, RI, May 1985
  5. 50 REM Random-entry matrices, single precision
  6. 60 REM Quartic modification of Newton's iteration
  7. 70 INPUT "Run number = ";RN
  8. 80 REM N data is in line 640, D1 data in line 660
  9. 90 READ N, D1:RANDOMIZE RN:OPTION BASE 1
  10. 100 REM Can optionally insert DEFDBL A,B,E,X instruction here
  11. 110 REM A is the matrix to invert, B is the approximate inverse of A, E is the          error matrix for the inverse, and X is a temporary storage matrix as            explained in the BYTE article.
  12. 120 DIM A(N,N),B(N,N),E(N,N),X(N,N)
  13. 130 REM Generate random-entry A-matrix. For real data, replace 140 by an INPUT          statement and input normalized data (each element of A is divided by the        largest element, L). Delete line 150. B and D1 must be multiplied by L          later.
  14. 140 FOR I=1 TO N:FOR J=1 TO N:A(I,J)=RND(1):NEXT :NEXT
  15. 150 FOR I=1 TO N:FOR J=1 TO N:IF RND(1)<.5 THEN A(I,J)=-A(I,J)
  16. 160 NEXT :NEXT
  17. 170 CLS:TIME$="00"   'Start clock
  18. 180 REM Eval. t by Pan-Reif eq 8 in BYTE article
  19. 190 FOR I=1 TO N:R0=0:S0=0:FOR J=1 TO N
  20. 200 R0=R0+ABS(A(I,J)):S0=S0+ABS(A(J,I)):NEXT J
  21. 210 X(1,I)=R0:E(1,I)=S0:NEXT I:T1=0:T2=0
  22. 220 FOR I=1 TO N:IF X(1,I)>T1 THEN T1=X(1,I)
  23. 230 IF E(1,I)>T2 THEN T2=E(1,I)
  24. 240 NEXT I:T=1/(T1*T2)
  25. 250 REM Eval. initial B-matrix
  26. 260 FOR I=1 TO N:FOR J=1 TO N:B(I,J)=T*A(J,I):NEXT :NEXT :H=0
  27. 270 PRINT "ITER."TAB(10);"E(1,1)"TAB(28);"E(1,N)"TAB(46);"B(1,1)"TAB(64);"B(1,N)"
  28. 280 REM Eval. error matrix E
  29. 290 FOR I=1 TO N:FOR K=1 TO N:Z=0:FOR J=1 TO N
  30. 300 Z=Z+B(I,J)*A(J,K):NEXT J:E(I,K)=-Z:NEXT K:NEXT I
  31. 310 FOR I=1 TO N:E(I,I)=1+E(I,I):NEXT I
  32. 320 BEEP:PRINT H;TAB(6);E(1,1);TAB(24);E(1,N);TAB(42);B(1,1);TAB(60);B(1,N)
  33. 330 IF H>50 THEN PRINT "STUCK!":BEEP:BEEP:GOTO 620
  34. 340 REM Test for escape from loop
  35. 350 FOR I=1 TO N:FOR J=1 TO N:IF ABS(E(I,J))>D1 THEN 370
  36. 360 NEXT :NEXT :GOTO 440
  37. 370 H=H+1  'Newton's iteration (quartic modification)
  38. 380 FOR I=1 TO N:X(1,I)=1+(1+(1+(1+E(I,I))*E(I,I))*E(I,I))*E(I,I):NEXT
  39. 390 FOR I=1 TO N:FOR J=1 TO N:E(I,J)=X(1,I)*E(I,J):NEXT :E(I,I)=1+E(I,I):NEXT
  40. 400 FOR I=1 TO N:FOR K=1 TO N:W=0:FOR J=1 TO N
  41. 410 W=W+E(I,J)*B(J,K):NEXT J:X(I,K)=W:NEXT K:NEXT I
  42. 420 FOR I=1 TO N:FOR J=1 TO N:B(I,J)=X(I,J):NEXT :NEXT
  43. 430 GOTO 280
  44. 440 REM output follows
  45. 450 T$=TIME$:BEEP:BEEP:BEEP:PRINT "Run number";RN,"Duration = ";T$
  46. 460 PRINT "No. iter.=";H,"Matrix dim.="N,"Max. error=";D1
  47. 470 INPUT "Another iteration (y,n)";C$:IF C$="Y" OR C$="y" THEN TIME$="00":GOTO 370
  48. 480 PRINT "The maximum error was calculated by multiplying A^-1*A."
  49. 490 PRINT "As a check on the true error in the approximated inverse matrix"
  50. 500 PRINT "the maximum error in A*A^-1 will now be determined."
  51. 510 ERASE X      'recycle the temporary storage matrix for errors in A*A^-1
  52. 520 DIM X(N,N)
  53. 530 FOR I=1 TO N:FOR J=1 TO N:FOR K=1 TO N
  54. 540 X(I,K)=X(I,K)+A(I,J)*B(J,K)
  55. 550 NEXT :NEXT :NEXT
  56. 560 REM You may print out the matrix and its inverse by placing a PRINT command         for A and B in the following loop.
  57. 570 FOR I=1 TO N:FOR J=1 TO N
  58. 580 TEST.VAL=ABS(X(I,J)):IF I=J THEN TEST.VAL=ABS(1-TEST.VAL)
  59. 590 IF TEST.VAL>MAX.ERROR THEN MAX.ERROR=TEST.VAL
  60. 600 NEXT :NEXT
  61. 610 PRINT "The max. error found this way is ";MAX.ERROR
  62. 620 END
  63. 630 REM Matrix dim. N =
  64. 640 DATA 10
  65. 650 REM Error criterion D1 =
  66. 660 DATA 3E-6
  67.